{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to use pythreejs to plot a superellipsoid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A superellipsoid is given by a parametric function and the equation is very similar to an ellipse equation. We only have different exponents which give us different shapes. For more informations: https://en.wikipedia.org/wiki/Superellipsoid.\n",
    "\n",
    "The idea of this example is to construct the mesh of the square $[0, 1]\\times[0,1]$ and to do a projection of these points on the superillipse which is the 2D shape and then to do a spherical product to have the 3D shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "n = 10 # number of discretisation points for the square in each direction \n",
    "x_box = np.concatenate((np.linspace(-1, 1., n), np.ones(n-2), np.linspace(1, -1., n), -np.ones(n-2)))\n",
    "y_box = np.concatenate((-np.ones(n-1), np.linspace(-1, 1., n), np.ones(n-2), np.linspace(1, -1., n-1, endpoint=False)))\n",
    "nx_box = x_box.size\n",
    "\n",
    "coords = np.empty((nx_box**2, 3), dtype=np.float32)\n",
    "\n",
    "def superellipse(rx, ry, m):\n",
    "    \"\"\"\n",
    "    superellipse formula with the projection of the unit square\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    rx : the radius in the x direction \n",
    "    ry : the radius in the y direction \n",
    "    m : the exponent of the superellipse\n",
    "    \n",
    "    Output\n",
    "    ------\n",
    "    the coordinates of the superellipse\n",
    "    \"\"\"\n",
    "    return x_box*rx*(1. - .5*np.abs(y_box)**(2./m))**(m/2.),  y_box*ry*(1. - .5*np.abs(x_box)**(2./m))**(m/2.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def superellipsoid(rx, ry, rz, m1, m2):\n",
    "    \"\"\"\n",
    "    superellipsoid formula with the spherical product of two superellipse\n",
    "    and update of the global coords array\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    rx : the radius in the x direction \n",
    "    ry : the radius in the y direction \n",
    "    rz : the radius in the z direction \n",
    "    m1 : the exponent of the first superellipse\n",
    "    m2 : the exponent of the second superellipse\n",
    "    \"\"\"    \n",
    "    gx, gy = superellipse(1, 1, m2)\n",
    "    hx, hy = superellipse(1, 1, m1)\n",
    "\n",
    "    \n",
    "    coords[:, 0] = rx*(gx[np.newaxis, :]*hx[:, np.newaxis]).flatten()\n",
    "    coords[:, 1] = ry*(gx[np.newaxis, :]*hy[:, np.newaxis]).flatten()\n",
    "    coords[:, 2] = rz*(gy[np.newaxis, :]*np.ones(hx.size)[:, np.newaxis]).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# superellipsoid parameters\n",
    "rx = ry = rz = 1.\n",
    "m1 = m2 = 1.\n",
    "\n",
    "superellipsoid(rx, ry, rz, m1, m2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct the triangulation by using the ConveHull function in scipy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.spatial as spatial\n",
    "\n",
    "cvx = spatial.ConvexHull(coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pythreejs import *\n",
    "from ipydatawidgets import ConstrainedNDArrayWidget\n",
    "from IPython.display import display\n",
    "\n",
    "view_width = 600\n",
    "view_height = 400\n",
    "\n",
    "coordinate_widget = ConstrainedNDArrayWidget(dtype='float32')(array=coords)\n",
    "\n",
    "surf_g = BufferGeometry(attributes=dict(\n",
    "    position=BufferAttribute(coordinate_widget),\n",
    "    index=BufferAttribute(cvx.simplices.ravel().astype(np.uint16)),\n",
    "))\n",
    "surf = Mesh(geometry=surf_g, material=MeshBasicMaterial(color='green', wireframe=True))\n",
    "scene = Scene(children=[surf, AmbientLight(color='#777777')])\n",
    "c = PerspectiveCamera(position=[2, 2, 3], up=[0, 0, 1],\n",
    "                      aspect=view_width / view_height,\n",
    "                      children=[DirectionalLight(color='white',\n",
    "                                                 position=[3, 5, 1],\n",
    "                                                 intensity=0.6)])\n",
    "renderer = Renderer(\n",
    "    camera=c, scene=scene, controls=[OrbitControls(controlling=c)], width=view_width, height=view_height)\n",
    "display(renderer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import FloatSlider, HBox, VBox\n",
    "\n",
    "m1_slider, m2_slider = (FloatSlider(description='m1', min=0.01, max=4.0, step=0.01, value=m1,\n",
    "                                            continuous_update=False, orientation='vertical'),\n",
    "                        FloatSlider(description='m2', min=0.01, max=4.0, step=0.01, value=m2,\n",
    "                                            continuous_update=False, orientation='vertical'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rx_slider, ry_slider, rz_slider = (FloatSlider(description='rx', min=0.01, max=10.0, step=0.01, value=rx, \n",
    "                                               continuous_update=False, orientation='horizontal'),\n",
    "                                   FloatSlider(description='ry', min=0.01, max=10.0, step=0.01, value=ry, \n",
    "                                               continuous_update=False, orientation='horizontal'),\n",
    "                                   FloatSlider(description='rz', min=0.01, max=10.0, step=0.01, value=rz, \n",
    "                                               continuous_update=False, orientation='horizontal'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update(change):\n",
    "    superellipsoid(rx_slider.value, ry_slider.value, rz_slider.value, \n",
    "                   m1_slider.value, m2_slider.value)\n",
    "    coordinate_widget.array = coords\n",
    "    coordinate_widget.notify_changed()\n",
    "    \n",
    "m1_slider.observe(update, names=['value'])\n",
    "m2_slider.observe(update, names=['value'])\n",
    "rx_slider.observe(update, names=['value'])\n",
    "ry_slider.observe(update, names=['value'])\n",
    "rz_slider.observe(update, names=['value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VBox([HBox([renderer, m1_slider, m2_slider]), rx_slider, ry_slider, rz_slider])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  },
  "widgets": {
   "state": {
    "5eb4554c95814fc6b83135c8b1a7a5ee": {
     "views": [
      {
       "cell_index": 8
      }
     ]
    },
    "b9b3e029e6554c0caee4a3cd691eb821": {
     "views": [
      {
       "cell_index": 12
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
